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1.0  MAJOR  ACCOMPLISHMENTS 

We  have  studied  the  temperature  dependence  of  the  frequency  shift  and  the  amplitude 
dependence  of  the  vibration  frequency  for  <100>  length  extension  modes  in  n-doped  silicon  (Si) 
microelectromechanical  (MEMS)  resonators  and  have  started  a  comparison  with  the 
experimental  data  from  Kenny’s  group.  The  theoretical  analysis  is  based  on  minimizing  the  free 
energy 


F=Fe+Fv;  Fe  =  \  dr  {zj  [n„(r)]  +  na  (r)  (D 


over  the  electron  densities  na(r )  in  different  valleys  a  =  1,2,3  for  a  given  strain  su(r)  with  the 
constraint  that  the  total  density  remains  constant,  8F/8na(r)  =  /i.  Parameters  are  the 
parameters  of  the  deformation  potential,  with  typical  values  ~10eV;  Fv  is  the  free  energy  of 
vibrations  in  the  absence  of  electron-phonon  coupling.  The  minimization  applies  where  the  rate 
of  intervalley  electron  transitions  largely  exceeds  the  typical  frequencies  of  the  studied  modes, 
which  is  the  case  in  the  experiment.  The  physical  idea  is  that  strain  shifts  the  electron  energy 
minima  with  respect  to  each  other,  making  the  populations  of  different  minima  different,  as 
illustrated  in  Figure  1.  In  turn,  the  difference  in  the  populations  leads  to  stress. 


Figure  1:  Strain  Shifts  Valleys,  Leading  to  a  Population  Difference 

The  minimized  electron  free  energy  takes  the  form  of  a  series  expansion  in  the  strain  tensor. 


min  Fe  =  \a2H  +  -  A 3iii  +  —A Aiiii  + 

2  6  4  24  4 


Tensors  An  are  proportional  to  the  electron  density  and  to  the  parameter  (swr  with 
n  =  2,3,4.  This  parameter  is  extremely  large  for  room  temperature  and  for  the  typical  doping 
densities  of  1019  cm'3 ,  reaching  the  value  of  order  of  103  »  1. 


At  the  beginning  of  the  project,  the  analysis  was  applied  to  <100>  extensional  modes,  where 
tensor  Sit  has  only  diagonal  components  along  the  axes  of  the  electron  valleys  and  the  coupling 
energy  — ik^ik  is  nonzero.  In  the  analysis  we  have  also  taken  into  account  the  internal 
nonlinearity  of  the  crystal  lattice,  which  leads  to  the  temperature  dependence  of  the  mode 
frequencies  even  in  the  absence  of  the  electron-phonon  coupling.  This  effect  is  additive  in  terms 


(2) 
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of  the  frequency  shift,  as  seen  from  the  expression  for  the  free  energy.  The  quartic  in  i  term  of 
internal  nonlinearity  (in  energy)  is  small  compared  to  the  electron-phonon  coupling  induced 
nonlinearity,  and  therefore  it  may  be  disregarded  in  calculating  the  nonlinear  frequency  shift. 

In  the  second  quarter,  we  extended  the  results  to  <1 10>  modes.  In  these  modes,  strain  tensor  has 
nonzero  off-diagonal  (shear)  components.  A  major  extension  of  the  project  is  that  we  have 
started  exploring  the  effective  two-phonon  coupling  that  results  from  the  phonon-induced 
anticrossing  of  the  bands.  In  silicon,  in  the  conventionally  used  effective  mass  approximation 
there  is  no  linear  coupling  between  shear  strain  and  conduction  electrons.  However,  shear  strain 
has  pronounced  effects  on  the  electronic  structure  beyond  this  approximation.  It  comes  primarily 
from  the  vicinity  of  the  degeneracy  (Dirac-type)  point  in  the  band  structure,  where  two  electron 
energy  bands  intersect.  This  is  point  X  on  the  boundary  of  the  Brillouin  zone  shown  in  Figure  2 
below. 


The  Hamiltonian  around  point  X  along  the  kz  axis  is  written  as: 

Hx  —  IA  H"  (TX(^A^ k\ky  "h  —  3^jcv)  "t”  ^4 


Here, 


A  A i kz  +  A2(kx  +  Ezz  +  ^2^xx  ^yy) 

As  a  result,  the  kz-  valley’s  bottom  shifts  quadratically  with  shear  strain: 

=  -%eiy 


(3) 

(4) 

(5) 


where  SE  is  the  splitting  of  the  energy  bands  at  the  positons  of  the  valleys.  A  similar  expression 
applies  to  other  valleys.  Effectively,  this  result  is  equivalent  to  a  two-phonon  coupling. 
Incorporating  this  effect  into  the  model  described  above  makes  it  possible  to  calculate  the  effects 
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of  the  shear  electron-phonon  coupling  on  the  linear  and  nonlinear  elasticity  for  <1 10>  modes. 
This  is  a  major  endeavor,  and  the  calculations  are  at  the  early  stage. 

In  the  third  quarter  we  have  performed  detailed  calculations  of  the  two-phonon  coupling  due  to 
shear  strain.  We  have  shown  how  it  affects  the  temperature  dependence  of  the  mode  frequencies 
and  the  mode  nonlinearity,  including  the  dependence  of  the  mode  frequency  on  the  vibration 
amplitude.  We  have  also  started  a  comparison  of  the  results  with  the  experiments  from  the 
Kenny  group  on  the  frequency  dependence  on  the  temperature  and  the  amplitude. 

In  the  fourth  quarter  we  have  summarized  the  results  in  a  paper,  which  is  included  as 
Appendix  A,  and  have  developed  a  code  that  is  publicly  available  at: 
http://www.pa.msu.edu/people/dykman/nonlinear  elasticity 

This  code  allows  calculating  the  temperature  dependent  corrections  to  the  frequency  due  to  the 
doping  and  the  nonlinear  frequency  shift  due  to  finite  vibration  amplitude  in  Si  resonators  for 
different  types  of  modes.  The  code  provides  both  the  analytical  expressions  and  evaluates  them 
numerically  using  the  parameters.  It  also  allows  finding  the  doping-induced  corrections  to  linear 
and  nonlinear  elasticity  in  germanium  (Ge). 
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2.0  MAJOR  FINDINGS 


•  Fully  analyzed  the  effect  of  the  degeneracy  of  the  conduction  bands  at  the  points  X  of  the 
Brillouin  zone  on  the  mode  nonlinearity  and  the  temperature  dependent  frequency  shift  due 
to  the  doping.  The  analysis  goes  beyond  the  conventional  deformation  potential 
approximation  and  is  necessary  to  describe  the  coupling  to  shear  strain  modes. 

•  Tested  the  theory  against  the  existing  in  the  literature  experimental  data  on  the  effect  of 
doping  on  the  speed  of  transverse  sound  and  the  hitherto  unexplained  effect  of  uniaxial  stress 
on  the  speed  of  sound  in  the  geometry  where  the  deformation  potential  does  not  lead  to  such 
an  effect. 

•  In  collaboration  with  Kenny’s  group,  obtained  new  experimental  data  on  the  doping-  induced 
temperature  dependence  of  the  frequency  and  on  the  nonlinearity  of  the  <1 10>  and  <100> 
Lame  modes. 

•  Described  the  experimental  data  on  the  temperature  and  amplitude  dependence  of  the 
frequency  of  the  Lame  modes  in  Si  microresonators.  The  results  describe  the  experimental 
observations  on  the  anomalous  softening  nonlinearity  in  the  <1 10>  modes,  which  behaves 
monotonically  with  doping  density,  as  opposed  to  the  hardening  nonlinearity  in  the  <100> 
modes  which  behaves  in  a  non-monotonic  way  with  doping  density. 

•  Developed  a  symbolic  and  numerical  codes  to  evaluate  the  effects  for  different  modes  in  Si 
and  Ge  microresonators. 

The  modes  used  in  the  simulations  of  Si  resonators  are  shown  in  Figure  3.  The  results  of  a 
comparison  of  the  theory  with  the  experimental  data  on  the  nonlinear  frequency  shift  of  the 
Lame  modes  in  Si  microresonators  are  shown  in  Figure  4.  The  comparison  is  performed  for  two 
donor  densities.  The  data  refer  to  two  types  of  modes.  One  mode  is  coupled  to  the  electron 
density  via  deformation  potential  (the  Lame  mode  in  a  plate  cut  along  the  <100>  axes).  The  other 
mode  refers  to  the  plate  cut  along  <1 10>  axes,  where  the  displacement  field  in  is  pure  shear  and 
the  coupling  to  the  electron  density  is  through  the  band  splitting  at  the  X  points  of  the  Brillouin 
zone. 


h  i 


Figure  3:  Two  Types  of  Vibrational  Modes  in  Microelectromechanical  Systems  (MEMS) 
used  to  calculate  the  Linear  and  Nonlinear  Frequency  Shift 

Left  panel:  Length  extension  mode  in  the  <11 0>  direction  in  a  thin  narrow  resonator.  Right 
panel:  Lame  mode  in  a  thin  plate  resonator  with  the  sides  along  <110>. 
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Figure  4:  Nonlinear  Frequency  Shift  Scaled  by  the  Squared  Vibration  Amplitude  dx  vs. 

Temperature  in  n-doped  Si  Lame  Resonators 

Solid  lines  show  the  present  theory.  Dots:  experimental  points  obtained  by  the  group  of 
T.  Kenny.  The  experimental  values  of  the  amplitude  are  obtained  assuming  the  gap  between  the 
resonator  and  the  electrodes  to  be  1.1 jum  in  the  absence  of  vibrations.  The  donor  density  in  the 
left  panel  is:  n  =  2.8  x  10I8cm3;  in  the  right  panel  n  =  5.9  x  1019  cm3. 

The  theory  has  no  adjustable  parameters  and  is  in  an  excellent  agreement  with  the  experiment. 

The  obtained  results  show  that  the  goal  of  the  project  has  been  reached.  We  have  full  quantitative 
theory  of  the  temperature-  and  density-dependent  change  of  the  vibration  frequency  of 
eigenmodes  in  resonators  and  of  their  dependence  on  the  vibration  amplitude. 
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3.0  THE  TEAM 


The  work  was  done  at  Michigan  State  University.  The  principal  investigators  were  Mark 
Dykman  and  Steven  Shaw.  The  graduate  student  was  Kirill  Moskovtsev.  The  work  was  done 
close  cooperation  with  the  team  of  Thomas  Kenny  at  Stanford  University. 
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APPENDIX  A:  STRONG  VIBRATION  NONLINEARITY  IN 
SEMICONDUCTOR-BASED  NANOMECHANICAL  SYSTEMS 

Kirill  Moskovtsev1  and  M.  I.  Dykman1 

department  of  Physics  and  Astronomy,  Michigan  State  University,  East  Lansing,  Michigan  48824,  USA 

(Dated:  November  28,  2016) 

We  study  the  effect  of  the  electron-phonon  coupling  on  vibrational  eigenmodes  of  nano-  and 
micro-mechanical  systems  made  of  semiconductors  with  equivalent  energy  valleys.  We  show 
that  the  coupling  can  lead  to  a  strong  mode  nonlinearity.  The  mechanism  is  the  lifting  of  the 
valley  degeneracy  by  the  strain.  The  redistribution  of  the  electrons  between  the  valleys  is 
controlled  by  a  large  ratio  of  the  electron-phonon  coupling  constant  to  the  electron  chemical 
potential  or  temperature.  We  find  the  quartic  in  the  strain  terms  in  the  electron  free  energy,  which 
determine  the  amplitude  dependence  of  the  mode  frequencies.  This  dependence  is  calculated  for 
silicon  micro-  systems.  It  is  significantly  different  for  different  modes  and  the  crystal  orientation, 
and  can  vary  nonmonotonously  with  the  electron  density  and  temperature. 

I.  Introduction 

The  electron-phonon  coupling  strongly  affects  vibrational  modes  of  nano-  and  micro-electro- 
mechanical  systems.  Much  interest  have  attracted  the  effects  of  this  coupling  related  to  the 
reduced  dimensionality  of  the  electron  system,  as  they  make  it  possible  to  reveal  interesting 
consequences  of  the  electron  correlations  at  the  nanoscale,  the  Coulomb  blockade  being  a  simple 
example,  cf.  [1-9]  and  references  therein. 

Much  less  attention  has  been  paid  to  the  consequences  of  the  electron-phonon  coupling,  which 
are  related  to  the  discreteness  of  the  vibrational  spectrum  of  a  nanosystem,  but  emerge  in  the 
absence  of  size  quantization  of  the  electron  motion.  One  of  such  consequences,  which  we  study 
in  this  paper,  is  the  coupling-induced  change  of  the  vibration  nonlinearity.  Strong  nonlinearity  is 
a  generic  feature  of  vibrations  in  small  systems  [10,  11],  Its  easily  accessible  manifestation  is  the 
dependence  of  the  mode  frequencies  on  the  vibration  amplitudes.  This  dependence  corresponds 
to  the  “self-action”  of  the  mode,  and  its  familiar  analog  in  bulk  crystals  are  acoustic  soli-  tons 
[12,  13];  however,  the  nonlinearity  required  for  observing  such  solitons  usually  is  sufficiently 
strong  only  for  high-frequency  phonons.  Also,  the  change  of  the  eigen  frequency  with  the  mode 
amplitude  is  of  interest  for  modes  with  a  discrete  frequency  spectrum,  such  as  standing  waves  in 
mesoscopic  systems,  but  not  for  propagating  waves  with  a  quasi-continuous  spectrum. 

Much  attention  have  been  recently  attracting  Si-based  nano-  and  micromechanical  systems,  see 
[14,  15]  and  references  therein.  In  such  systems  there  was  observed  an  unexpectedly  large 
change  of  the  amplitude  dependence  of  the  vibration  frequency  with  the  varying  electron  density 
[16,  17],  When  the  doping  level  was  increased  from  2.8  x  1018  cm-3  to  5.9  x  1019  cm-3,  the 
nonlinearity  parameter  increased  by  more  than  an  order  of  magnitude.  Moreover,  the  nonlinearity 
change  was  different  for  the  vibrational  modes  with  different  spatial  structure.  In  this  paper  we 
develop  a  theory  of  the  nonlinearity  of  vibrational  modes  in  semiconductor  nano-  and  micro¬ 
mechanical  systems  with  high  electron  density.  We  show  that  the  electron-phonon  coupling  can 
lead  to  a  strong  self-action  of  the  vibrational  modes,  which  in  turn  significantly  modifies  the 
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amplitude  dependence  of  the  mode  frequencies.  We  find  the  dependence  of  the  effect  on  the 
electron  density  and  temperature. 

For  bulk  semiconductors,  the  effect  of  the  electron-  phonon  coupling  on  the  elastic  properties, 
including  the  three-phonon  coupling,  was  first  analyzed  by  Keyes  [18].  The  analysis  referred  to 
n-Ge  and  was  based  on  the  deformation  potential  approximation.  The  idea  was  that  de-  formation 
lifts  the  degeneracy  of  the  equivalent  electron  valleys,  which  leads  to  a  redistribution  of  the 
electrons  over  the  valleys.  In  turn,  such  redistribution  changes  the  speed  of  sound  depending  on 
the  direction  and  polarization  of  the  sound  waves  and  also  affects  the  sound  speed  in  the 
presence  of  uniaxial  stress.  This  theory  was  ex-  tended  to  silicon  and  the  corresponding 
measurements  were  done  by  Hall  [19].  However,  Hall  also  observed  the  change  of  the  speed  of 
transverse  sound  waves  and  the  effect  of  stress  on  sound  propagation  in  the  geometries,  where 
these  effects  are  due  to  shear  deformation  and  do  not  arise  in  the  deformation  potential  model.  A 
theory  of  the  change  of  the  linear  shear  elastic  constant  in  silicon  due  to  the  intervalley 
redistribution  of  the  electrons  was  developed  by  Cerdeira  and  Cardona  [20]. 

As  we  show,  in  mesoscopic  systems  the  strain-induced  redistribution  of  the  electrons  over  the 
valleys  of  the  conduction  band  leads  to  the  previously  unexplored  strong  fourth-order 
nonlinearity  of  the  vibrational  modes.  This  nonlinearity  gives  a  major  contribution  to  the 
amplitude  dependence  of  the  vibration  frequency.  The  redistribution  also  leads  to  a  temperature 
dependence  of  the  frequencies.  The  magnitudes  of  the  effects  sensitively  de-  pends  on  the  mode 
structure.  We  describe  them  for  several  types  of  modes,  including  those  studied  in  the 
experiment  [16,  17]  and  qualitatively  compare  the  results  with  the  observations.  The  theoretical 
results  refer  to  both  degenerate  and  nondegenerate  electron  systems.  Specific  calculations  are 
done  for  silicon  resonators. 

In  Sec.  II  we  give,  for  completeness,  the  expressions  for  the  mode  normalization  and  the 
amplitude-dependent  frequency  shift  of  coupled  nonlinear  modes  in  a  nano-  or  micro-system.  In 
Sec.  Ill  and  Appendix  A  we  provide  expressions  for  the  electron-phonon  coupling  induced 
change  of  the  elasticity  parameters,  including  the  parameters  of  quartic  nonlinearity.  In  Sec.  IV 
we  discuss  the  asymptotic  behavior  of  the  parameters  of  quartic  nonlinearity  for  low  and  high 
electron  density  and  give  their  explicit  form  for  silicon.  In  Sec.  V  we  calculate  the  nonlinear 
frequency  shift  for  several  frequently  used  vibrational  modes  in  single-crystal  silicon  systems 
and  show  the  dependence  of  this  shift  on  the  electron  density  and  temperature.  The  explicit 
analytical  expressions  are  given  in  Appendices  C  and  D.  Sec.  VI  contains  concluding  remarks. 

II.  Nonlinear  Frequency  Shift  of  Low-Frequency  Eigenmodes 

Of  primary  interest  for  nano-  and  micro-mechanical  systems  are  comparatively  low-frequency 
modes  with  wavelength  on  the  order  of  the  maximal  size  of  the  system.  Examples  are  provided 
by  long-wavelength  flexural  modes  of  nanotubes,  nanobeams,  and  nano/micro-membranes,  or 
acoustic-type  modes  in  microplates  or  beams.  These  modes  are  easy  to  excite  and  detect.  We  will 
enumerate  them  by  index  v.  Their  dynamics  is  de-  scribed  by  the  elasticity  theory  [21].  The 
spatial  structure  of  the  displacement  field  of  a  mode  u(v)(r)  in  the  harmonic  approximation  is 
determined  by  the  boundary  conditions.  We  will  choose  u(v)(r)  dimensionless,  so  that  in  our 
finite-size  system 
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Here,  V  is  the  volume  of  the  system.  We  assumed  that  the  mode  eigenfrequencies  cov  are 
nondegenerate;  including  degenerate  modes  is  straightforward.  For  simplicity,  we  also  assumed 
that  the  system  is  spatially  uniform;  an  extension  to  spatially  nonuniform  systems  is 
straightforward  as  well. 

We  emphasize  the  distinction  of  the  normalization  (1)  from  the  conventional  normalization  for 
bulk  crystals,  where  v  corresponds  to  the  wave  vector  and  the  branch  number,  and  the 
normalization  integral  is  independent  of  the  volume.  The  normalization  (1)  is  convenient  for  the 
analysis  of  low-frequency  modes  with  the  discrete  spectrum  characteristic  of  mesoscopic 
systems.  Such  modes  are  standing  waves,  and  therefore  vectors  u(v)  can  be  chosen  real. 

The  low-frequency  part  of  the  displacement  can  be  written  as 


Functions  Qv{t )  give  the  mode  amplitudes.  In  the  harmonic  approximation  the  dynamics  of  the 
standing  waves  is  described  by  the  Hamiltonian 


where  Pv  is  the  momentum  of  mode  v  and  M  is  the  mass  of  the  system. 

The  anharmonicity  of  the  crystal  leads  to  mode-mode  coupling.  Within  the  elasticity  theory  this 
coupling  is  de-  scribed  by  the  terms  in  the  Hamiltonian,  which  are  cubic  and  quartic  in  the  strain 
tensor.  We  will  not  consider  higher-order  terms,  which  are  small  for  the  mode  amplitudes  of 
interest.  From  the  expansion  (2),  we  obtain  the  nonlinear  part  of  the  Hamiltonian  in  the  form 


Equation  (4)  is  essentially  an  expansion  in  the  ratio  of  the  mode  amplitudes  to  their  characteristic 
wavelength,  which  is  of  the  order  of  the  appropriate  linear  dimension  of  the  system.  This  is  why 
mesoscopic  systems  are  of  particular  interest,  as  here  vibrations  of  low-frequency  eigenmodes 
become  nonlinear  for  already  small  vibration  amplitudes. 

A  familiar  consequence  of  nonlinearity  in  nano-  and  micromechanical  systems  is  the  dependence 
of  the  vibration  frequency  of  a  mode  on  its  own  amplitude  and  on  the  amplitudes  of  other  modes. 
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see  Ref.  1 1  for  a  review.  In  particular,  the  change  deov  of  the  mode  frequency  due  to  the 
vibrations  of  the  mode  itself,  Qv  (t)  =  Av  cos  covt,  is  [22,  23] 


SbJv  S3 ■ 


3-fi,  y-  -S.^} 

8\l^f 


(*) 


where  fv  —  yww  and  we  kept  the  terms  of  the  first  order  in  y  and  the  second  order  in  /?. 

The  nonlinear  mode  coupling  (4)  leads  also  to  the  frequency  shift  due  to  thermal  vibrations  of  the 
modes.  The  dominating  contribution  to  this  shift  for  low-frequency  modes  comes  from  their 
coupling  to  modes  with  frequencies  ~  kuT/h.  which  have  a  much  higher  density  of  states.  This 
shift  is  described  by  an  expression  that  is  similar  to  Eq.  (5)  wish  replaced  by  A?jt  ~  ksT /Muj 
and  placed  under  the  sum  over  v’,  in  the  classical  limit. 

III.  The  Nonlinearity  Due  to  the  Electron-Phonon  Coupling 

We  will  consider  the  vibration  nonlinearity  due  to  the  electron-phonon  coupling  in  multi-valley 
semiconductors  with  cubic  symmetry,  silicon  and  germanium  being  the  best  known  examples.  In 
such  semiconductors,  the  energy  valleys  of  the  conduction  band  are  located  at  high-  symmetry 
axes  of  the  Brillouin  zone.  Strain  lifts  the  symmetry  and  thus  the  degeneracy  of  the  valleys. 

The  simplest  mechanism  of  the  electron-phonon  coupling  is  the  deformation  potential.  Here,  the 
energy  shift  8Ea  of  valley  a  is  determined  by  the  deformation  potential  parameters  Eu  and  Ed  of 
the  coupling  to  a  uniaxial  strain  along  the  symmetry  axis  of  the  valley  and  to  dilatation, 
respectively.  In  terms  of  the  stain  tensor  sy  we  have  _.£•?>£.  where 

®  e<Q),  with  being  the  unit  vector  along  the  symmetry  axis  of  the 
valley.  We  use  the  hat  symbol  to  indicate  tensors  and  symbol  “®”  to  indicate  tensor  products. 
The  analysis  below  is  not  limited  to  the  deformation  potential  approximation.  An  important 
extension  will  be  discussed  using  silicon  as  an  example. 

We  assume  that  the  strain  varies  in  time  and  space  slowly  compared  to  the  reciprocal  rate  of 
intervalley  electron  scattering  and  the  intervalley  scattering  length,  respectively.  Then  the 
electron  system  follows  the  strain  adiabatically.  The  electron  density  nia](r)  in  valley  a  is 
decreased  or  increased  depending  on  whether  the  bottom  of  the  valley  goes  up  or  down.  In  the 
single-electron  approximation  and  for  the  deformation  potential  coupling,  the  electron  free 
energy  density  for  a  given  strain  is  pe  _  4.  n*“>(r)s£^s  i(r)}  where  fe[n(r)]  is 

the  free  energy  density  for  electrons  with  density  n{ r)  in  a  valley  in  the  absence  of  coupling  to 
phonons. 

The  electro-neutrality  requires  that  the  total  electron  density  summed  over  the  valleys  be 
constant.  The  free  energy  density  Fe  has  to  be  minimized  over  n<a>(r)  to  meet  this  constraint. 

This  gives  the  change  of  the  electron  chemical  potential  8ju  due  to  strain  s'.  The  resulting 
increment  of  the  electron  free  energy  density  has  the  form  of  a  series  expansion  in  the  strain 
tensor, 
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Here  Ai;  A2 ,  A3,  and  are  tensors  of  ranks  2,  4,  6,  and  8,  respectively.  They  are 
contracted  with  the  tensor  products  of  the  strain  tensor  s'.  Respectively,  At  are  the  electronic 
contributions  to  the  linear  (for  k  =  2)  and  nonlinear  (for  k  >  2)  elasticity  parameters  of  the 
crystal.  These  contributions  are  isothermal,  but  since  the  change  of  the  mode  frequencies  from 
the  electron-phonon  coupling  is  small  and  the  nonlinearity  is  also  small,  the  difference  with  the 
adiabatic  expressions  can  be  disregarded. 

To  the  third  order  in  if  the  expression  for  SFe  in  terms  of  the  shift  of  the  valleys  was  found  by 
Keyes  [18]  in  the  analysis  of  sound  wave  propagation.  However,  to  find  the  parameters  of  the 
quartic  nonlinearity  of  resonant  modes  in  small  systems,  which  is  of  primary  interest  to  us,  we 
also  need  to  keep  quartic  terms  in  Eq.  (6). 

As  seen  from  the  explicit  form  of  the  parameters  of  the  expansion  (6)  given  in  Appendix  A, 

A*  ix  k  =  1,2,. . .),  where  juo is  the  electron  chemical  potential  in 

the  absence  of  strain;  it  is  determined  by  the  total  (summed  over  the  valleys)  electron  density  n. 
Of  central  importance  for  the  analysis  is  that  parameter  aj  max(/./o,  kuT )  ~  103  for  electron 
densities  n  ~  1019  cm-3  and  room  temperatures,  i.e. 


Su/  nmx(^o,  kg/T)  !S>  1.  (7) 


As  a  consequence,  the  coefficients  at  the  nonlinear  in  e*  terms  in  Eq.  (6)  quickly  increase  with 
the  increasing  order  of  the  nonlinearity  [the  overall  series  (6)  is  converging  fast  because  of  the 
smallness  of  the  strain  tensor]. 

The  increase  of  A*,  with  k  allows  us  to  keep  in  if  only  the  terms  linear  in  the  lattice 
displacement,  i.e.,  to  set  ?;;/  =  ( \/2)(dui/dxj  +  du/dxi),  where  where  m  and  x,  are  the  components  of 
the  displacement  and  the  coordinates,  respectively.  Indeed,  in  this  case  a  kth  term  of  the  series 
(6)  is  of  order  k  in  the  displacement.  If  we  included  the  quadratic  in  dui/dxj  term  into  one  of  the  s' 
tensors  in  the  kth  term,  this  term  would  become  of  order  k  +  1  in  the  displacement.  However,  for 
linear  if  the  (k  +  l)th  term  in  the  series  (6)  is  also  of  the  (k  +  l)th  order  in  the  displacement,  but 
is  larger  by  factor  E«/max(jUo,  ksT). 


For  linear  if,  the  total  strain  is  a  sum  of  partial  contributions  of  strain  from  individual  modes. 

For  mode  v,  such  partial  contribution  is  expressed  in  terms  of  the  scaled  displacement  u(v)(r)  [see 


Eq.  (2)]  as  A  =  Qv/V),  where  e4jJ(r)  =  ^[du^i^/dxj +duf'>(r)/dxl] 

contrast  to  the  dimensionless  strain  tensor  if,  tensor  s"<v>  has  dimension  [length]"1 


We  note  that,  in 


From  Eq.  (6)  we  find  the  electronic  contributions  to  the  nonlinearity  parameters 
in  Hamiltonian  (4), 
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dr  A3 


dr  A4 
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£(l,|)  ®  (g,  el"’)  ®  £(*'*),  (8) 


where  s'<v>  =  s'<v>(f):  tensors  A/,,  are  independent  of  r.  Similarly,  the  electronic  contribution  to 
the  eigenfrequency  is 


jW  = 

V 


—  f 

2Mu„  J 


dr  A,  •  £(">  0  eM 


Generally  the  term  oc  A 2  leads  to  mode  mixing;  however,  if  the  mode  frequencies  are 
nondegenerate,  this  mixing  is  weak  and  can  be  disregarded,  to  the  leading  order  in  the  electron- 
phonon  coupling.  One  can  see  that  the  effect  of  the  static  stress  oc  Ai  can  be  disregarded  as 
well. 


The  jfequency  change  (9)  depends  on  temperature  because  of  the  temperature  dependence  of 
oc  A2.  The  nonlinearity  (8)  also  leads  to  a  temperature  dependence  of  the  mode 
eigenfrequency.  Together  they  modify  the  temperature  dependence  of  the  mode  eigenfrequencies 
compared  to  that  of  undoped  crystals.  This  modification  often  weakens  the  temperature 
dependence  of  the  eigenfrequencies,  which  proves  very  important  for  applications  of  micro¬ 
mechanical  systems  in  devices  that  work  in  a  broad  temperature  range  [24]. 

Equations  (6)-(9)  are  generic  and  apply  beyond  the  deformation  potential  approximation.  This  is 
of  particular  importance  for  silicon.  Here,  the  electron  band  valleys  lie  on  the  (lOO)-axes  close  to 
the  X-points  on  the  zone  boundaries  where  two  electron  energy  bands  cross.  Lattice  strain  can 
lead  to  a  band  splitting  at  X-points  and  a  shift  of  the  valleys  [25,  26],  Importantly,  this  shift 
results  from  a  shear  strain,  which  does  not  lead  to  a  linear  in  the  strain  shift  in  the  deformation 
potential  approximation.  The  valley  shift  is  quadratic  in  s'  in  this  case,  as  explained  in  Appendix 
A,  which  corresponds  to  an  effectively  two-phonon  coupling.  The  coupling  parameter  Ssh  is 
quadratic  in  the  strain-induced  band  splitting,  see  Eq.  (A2).  It  is  large,  much  larger  than  the 
constant  H„.  Therefore  the  arguments  given  below  Eq.  (7)  apply  in  this  case  as  well.  For  purely 
shear  strain  in  silicon,  terms  of  odd  order  in  s'  in  SFe,  Eq.  (6),  vanish. 

IV.  Explicit  Form  of  the  Tensors  of  Nonlinear  Elasticity 

Tensors  Ancan  be  obtained  by  minimizing  the  free  energy  density  of  the  electron  system  for  a 
given  strain  and  expanding  the  result  in  a  series  in  s'.  A  general  procedure  that  allows  one  to  find 
the  components  An  for  n  <  4  is  described  in  Appendix  A.  Using  the  symmetry  arguments,  the 
elasticity  tensors  are  conveniently  written  in  the  contracted  (Voigt)  notation  where  the  symmetric 
strain  tensor  is  associated  with  a  six-component  vector.  Then  the  nonlinear  elasticity  tensors  Ag 
and  X  become  tensors  of  rank  three  and  four  in  the  corresponding  vector  space.  We  use 
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notation  Sc"  for  tensors  A  in  these  notations  to  emphasize  that  we  are  calculating  corrections  to 
the  nonlinear  elasticity  tensors  due  to  the  electron-phonon  coupling. 

The  explicit  expressions  for  the  nonlinear  elasticity  tensors  Sc"  are  given  in  Table  I.  They  refer  to 
silicon  and  include  the  contributions  that  come  from  both  the  deformation  potential  coupling  and 
from  the  splitting  of  the  electron  bands  due  to  shear  strain.  In  the  deformation  potential 
approximation,  the  components  of  the  third-rank  tensor  Sc",  which  determine  the  cubic  in  the 
strain  terms  in  the  free  energy,  were  found  earlier  [19].  Therefore  we  give  only  the  components 
that  contain  a  contribution  from  shear  strain. 

The  fourth-rank  tensor  Sc"  determines  the  quartic  in  the  strain  terms  in  the  free  energy  and  has 
not  been  dis-  cussed  before,  to  the  best  of  our  knowledge.  We  give  all  independent  components 
of  this  tensor.  It  is  expressed  in  terms  of  the  derivative  of  the  electron  density  n  over  the 
chemical  potential  in  the  absence  of  strain  ju o,  which  is  a  familiar  thermodynamic  characteristic. 

It  is  intuitively  clear  that  the  considered  effect  of  the  change  of  the  electron  density  in  different 
valleys  in  response  to  strain  should  be  related  to  the  derivative  dn/d^o.  Interestingly  because  we 
consider  nonlinear  response  to  strain,  the  expressions  in  Table  I  contain  also  higher-order 
derivatives  of  n  over  juo.  As  we  will  see,  this  leads  to  a  nontrivial  behavior  of  the  nonlinear 
frequency  shift  with  varying  temperature  and  density.  The  considered  mechanism  of  the  strain- 
induced  inter-valley  electron  redistribution  does  not  contribute  to  the  components  0123  and  a 456, 
therefore  Sewn  =  Sc  1456  =  0. 

A.  Nonlinear  Elasticity  in  the  Limiting  Cases 

The  expressions  for  Sc"  simplify  in  the  case  of  low  doping  (or  high  temperature),  where  the 
electron  gas  is  strongly  nondegenerate,  and  in  the  opposite  case  of  a  strongly  degenerate  electron 
gas.  For  a  nondegenerate  gas,  where  the  chemical  potential  in  the  absence  of  strain  is 
f jo  <  0,  IjC/ol  »  ksT,  we  have  in  Table  I  Fm(x)  =  with  x  =  /Ji/ksT.  The  /ro-dependent 

factors  exp(/./o/  knT)  in  F\n  and  its  derivatives  cancel  each  other  in  the  expressions  for  Sc"  and 
drop  out  from  these  expressions.  The  dependence  of  Sc"  on  density  is  then  just  linear,  <5£°c  n. 
Parameters  Ci,...,4  in  Table  I  depend  only  on  temperature,  C 1  oc T~l,  C2  oc 7^3,  C3  ocT^2  and 

c4  ocr'. 

The  decrease  of  the  nonlinear  elasticity  parameters  with  increasing  temperature  in  a 
nondegenerate  electron  gas  is  easy  to  understand.  The  effect  we  consider  is  determined  by  the 
competition  between  the  energetically  favorable  unequal  population  of  the  electron  energy 
valleys  in  a  strained  crystal  and  the  entropically  more  favorable  equal  valley  population.  With 
increasing  temperature  the  entropic  factor  becomes  stronger,  leading  to  a  smaller  population 
difference  and  thus  smaller  effect  of  the  electron  system  on  the  vibrations. 

For  strong  doping,  where  /ro/ZcsT  »  1,  we  have  /jo  °°  n2/3,  and  then 

Fi/cj(x)  ~  fa?3'2  with  x  =  fiQ/kijT  .  Therefore  parameters  Ci,...,4in  Table  I  become 
temperature  independent,  with  nC\  ocnm,  11C2  ocn~l,  nC3  ocn~m,  and  nQ  ocnU3. 
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The  results  on  the  asymptotic  behavior  of  the  corrections  to  nonlinear  elasticity  are  not  limited  to 
silicon.  Since  parameters  C  1,23,4  are  given  by  the  coefficients  in  the  general  expansion  of  the  free 
energy  in  strain,  (Al),  these  results  can  be  applied  to  the  nonlinear  elasticity  induced  by  the 
electron-phonon  coupling  in  other  multi-  valley  semiconductors.  To  illustrate  this  point,  in 
Appendix  B  we  give  Sc"  tensor  in  germanium. 

The  difference  between  the  asymptotic  behavior  of  the  tensors  Sc"  in  the  limits  of  nondegenerate 
and  strongly  degenerate  electron  gas  can  lead  to  a  peculiar  density  and  temperature  dependence 
of  the  nonlinear  frequency  shift  of  the  vibrational  modes.  It  comes  from  the  coefficients  C 
containing  higher-order  derivatives  of  n  with  respect  to  /no-  In  the  transition  region  /./o  ~  kuT , 
thinking  of  the  competition  between  the  entropic  and  energetic  factors  does  not  provide  a  simple 
insight  into  the  behavior  of  Sc",  as  both  the  energy  and  the  entropy  are  complicated  functions  of 
density  and  temperature. 

TABLE  I.  The  change  of  the  components  of  the  nonlinear  elasticity  tensors  due  to  the  strain- 
induced  electron  redistribution  between  equivalent  energy  valleys  in  doped  silicon.  The 
coordinate  axes  are  chosen  along  the  (100)  axes.  Parameter  ESh  characterizes  the  effectively  two- 
phonon  coupling  to  shear  strain.  This  parameter  as  well  as  function  Fmix)  are  defined  in 
Appendix  A;  x  =  ho/JcbT  and  n  is  the  electron  density. 


<^144  =  —  2<5cj55 


<fcim  =  — 2<5cui2  =  2£cii22 
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;  -u  —  sh  Li 
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Ci  =  F[/2/Fx/2ki}T  =  dlnn/ef/io 

C2  =  (kBT)-aF[j2  [d2(l/F[,2)/dx2]  / F\ ,2 
=  (dn/dfio)2  [<P(dno/dn)/dfil]  fn 
C3  =  F[’/2/Fi,2{kBT)2  =  n~ 1  d2  n /dfi2 
Cj  =  F'Xj2/F\ /2knT  =  d\nn/dfio 


V.  Doping-Induced  Nonlinearity  of  Simple  Vibrational  Modes 

The  nonlinear  elasticity  tensors  in  Table  I  give  the  doping-induced  contributions  to  the 
nonlinearity  parameters  of  the  eigenmodes  of  micro-  and  nano-mechanical  systems.  These 
contributions  are  described  by  Eq.  (8).  As  mentioned  before,  an  important  characteristic  of  the 
mode  nonlinearity  is  the  dependence  of  the  mode  frequency  on  the  vibration  amplitude.  To  the 
leading  order,  it  is  given  by  Eq.  (5).  This  dependence  has  a  contribution  from  the  nonlinearity  of 
an  undoped  crystal,  which  is  quadratic  in  the  parameters  of  the  cubic  nonlinearity;  for  example, 
if  the  latter  is  described  by  the  Griineisen  constant,  the  corresponding  contribution  is  quadratic  in 
this  constant.  It  is  typically  small.  There  is  also  a  contribution  from  the  quartic  nonlinearity;  the 
parameters  of  such  nonlinearity  are  not  known  in  undoped  crystals  and  are  not  expected  to  be 
large.  Respectively,  the  amplitude  dependence  of  the  vibration  frequency  for  low-frequency 
modes  in  weakly  doped  single-crystal  micro-mechanical  systems  is  relatively  weak  [17]. 

A  feature  of  the  doping-induced  nonlinearity  described  by  Table  I  is  that  the  quartic  in  the  strain 
term  in  the  free  energy  has  a  large  coefficient  compared  to  the  cubic  term,  cf.  Eq.  (7)  and  the 
discussion  below  this  equation.  Therefore,  in  Eq.  (5)  for  the  amplitude  dependence  of  vibration 
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frequency  one  can  keep  only  the  Duffing  non-linearity  constant  yv.  The  contribution  from  the 
cubic  nonlinearity  terms  ocf i2vvv7  can  be  disregarded.  For  a  mode  v,  the  doping-induced 
contribution  to  yv  is  equal  to  7^1/1/  in  Eq.  (8). 

To  find  the  dependence  of  the  mode  frequency  on  the  vibration  amplitude  we  go  through  the 
following  steps.  First,  we  find  the  normal  modes  of  interest  for  the  given  geometry  of  the  system, 
with  account  taken  of  the  boundary  conditions,  and  normalize  the  displacements  u(v)(r)  as 
indicated  in  Eq.  (1).  We  use  u(v)(r)  to  find  the  strain  tensor  e"(v)(r).  The  result  is  substituted  into 
Eq.  (8)  and  is  convoluted  with  tensor  A4,  giving  the  value  of  yv  which  is  then  used  in  Eq.  (5)  to 
find  the  frequency  dependence  on  the  vibration  amplitude  Scov.  Of  particular  interest  is  the 
relative  frequency  shift  S(Ov/cov.  To  find  this  shift  to  the  leading  order,  one  can  disregard 
nonlinearity  when  calculating  the  eigenfrequency  cov.  Then,  from  Eq.  (5), 


_  _ 37 VA% _ 

8  f  dr  flip  -eM®  eW 


(10) 


where  AT  '  is  the  full  tensor  of  linear  elasticity,  which  includes  the  major  term  of  the  linear 
elasticity  of  the  undoped  crystal  and  the  doping-induced  correction  A$ . 


An  important  feature  of  the  relative  shift  Scov/cov  is  its  scaling  with  the  size  of  the  system.  The 
vibration  amplitude  Av  in  Eq.  (10)  can  be  scaled  by  the  lateral  dimension  L,  for  example  the 
length  of  a  nanobeam  or  a  nanowire  for  an  extension  mode,  or  the  size  of  the  square  for  a  Lame 
mode,  or  the  diameter  of  a  disk  for  a  breathing  mode  in  a  disk.  Respectively,  we  write  Av  =  rjvL. 
Then,  if  one  takes  into  account  the  explicit  form  (8)  of  the  parameter  7 v  =  jiVvu  ,  one  finds 
from  Eq.  (10)  that  the  ratio  6av/ (r]v2  av)  is  independent  of  the  system  size  for  the 
aforementioned  modes.  In  this  estimate  we  used  that  the  tensors  A  are  material  parameters  and 
are  independent  of  the  geometry.  We  also  used  that  modes  of  interest  have  typical  wavelength 
~L,  therefore  £tu}  scales  as  L"1. 


Most  of  the  experiments  in  nano-  and  micromechanics  are  done  with  nanobeams,  nanowires, 
membranes,  or  thin  plates.  In  such  systems  the  thickness  is  much  smaller  than  the  length  or,  in 
the  case  of  membranes  or  plates,  the  lateral  dimensions.  Then,  from  the  boundary  condition  of 
the  absence  of  tangential  stress  on  free  surfaces  [21],  it  follows  that  the  strain  tensor  s'  weakly 
depends  on  the  coordinate  normal  to  the  surface.  This  simplifies  the  denominator  in  Eq.  (10), 
making  it  proportional  to  the  thickness.  Similarly,  from  Eq.  (8)  yv  is  also  proportional  to  the 
thickness,  and  the  thickness  drops  out  of  Eq.  (10).  The  explicit  expressions  for  MwJ  and  yv  that 
deter-  mine  the  denominator  and  the  numerator  in  Eq.  (10),  respectively,  are  given  in  Appendices 
C  and  D  for  Lame  and  extension  modes.  These  expressions  are  cumbersome,  and  it  is 
convenient  to  use  symbolic  programming  to  obtain  them.  [27] 
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A.  Temperature  and  Electron  Density  Dependence  of  the  Scaled  Nonlinear  Frequency 
Shift 

The  scaled  ratio  Suv / that  characterizes  the  relative  nonlinear  frequency  shift  is  shown 
in  Figure  1  for  several  modes  that  are  often  used  in  single-crystal  silicon  MEMS.  This  ratio 
depends  on  the  type  of  the  mode  and  the  crystal  orientation.  Figure  1  refers  to  high-symmetry 
crystal  orientations,  in  which  case  the  modes  have  a  comparatively  simple  spatial  structure  and 
the  surfaces  can  be  made  smooth.  We  used  the  values  SM  =  8.8  eV  [28],  ESh  =  300  eV,  the 
effective  mass  for  density  of  states  meff  =  0.32 me  [29],  and  the  temperature-dependent  linear 
elasticity  parameters  given  in  Ref.  [30], 

Figure  1  shows  that  the  electron-redistribution  induced  nonlinearity  of  vibrational  modes  is  very 
strong.  For  the  ratio  of  the  vibration  amplitude  to  the  system  size  //  ~  10"4  and  the  mode 
eigenfrequency  oh/2n  ~  10  MHz,  the  frequency  change  can  be  as  large  as  Scov/ln  ~  0. 1  kHz. 
This  explains,  qualitatively,  the  observations  [17].  A  quantitative  comparison  with  the 
experiment  [17]  is  complicated,  as  the  observations  refer  to  different  samples.  Our  preliminary 
results  show  an  excellent  quantitative  agreement  with  the  data  obtained  for  the  same  sample  at 
different  temperatures  and  for  different  types  of  modes  [31]. 


FIGURE  1 .  Relative  change  da>v  /cov  of  the  vibration  frequency  of  a  mode  with  the  vibration 
amplitude  rjv  scaled  by  the  relevant  size  of  the  system,  cf.  Eq.  (10).  The  results  refer  to  single 
crystal  silicon  resonators.  Curves  1  and  2  refer  to  the  first  Lame  mode  in  square  plates  cut  in 
(100)  and  (110)  directions,  respectively.  In  this  case,  the  size  of  the  resonator  is  the  length  of  the 
side  of  the  square.  Curves  3  and  4  refer  to  the  first  extension  mode  in  beams  cut  in  (100)  and 
(110)  directions,  respectively.  In  this  case,  the  size  of  the  resonator  is  the  length  of  the  beam. 

The  nonlinear  frequency  shift  displays  several  characteristic  features,  as  seen  from  Figure  1 .  One 
of  them  is  the  strong  dependence  of  the  shift  on  the  type  of  the  mode  and  the  crystal  orientation. 
For  both  the  Lame  and  the  extension  mode,  the  shift  is  much  stronger  for  crystals  cut  out  in 
(100)  direction  than  in  (110)  direction.  This  is  a  consequence  of  the  electron  energy  valleys  lying 
along  the  (100)  axes,  making  the  system  more  “responsive”  to  the  lattice  displacement  along 
these  axes.  Interestingly,  in  the  both  configurations  the  shifts  for  the  Lame  modes  are  larger  than 
for  the  extension  modes. 
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A  somewhat  unexpected  feature  is  the  nonmonotonic  dependence  of  the  nonlinear  frequency 
shift  on  the  electron  density  and  temperature.  The  nonmonotoncity  occurs  in  the  range  where  the 
electron  system  is  close  to  degeneracy,  /uo/IcbT  ~  1,  and  it  strongly  depends  on  the  crystal 
orientation.  It  is  much  stronger  for  crystals  cut  in  (100)  then  (110)  directions.  For  a  crystal  cut  in 
(110)  direction,  both  the  density  and  temperature  dependence  of  the  shift  are  monotonic  in  the 
case  of  the  Lame  mode,  whereas  for  the  extension  mode  the  nonmonotonicity  is  weak. 

The  nonmonotonicity  of  the  frequency  shift  stems  from  the  behavior  of  the  parameters  nC2,3,4  in 
the  range  /no  ~  kuT .  As  seen  from  Table  I,  parameter  nCi  exponentially  increases  with  the 
increasing  /Jo/kuT  for  negative  /Jo/kBT,  but  for  large  positive  /Jo/knT  it  falls  off  /ro//:«T3/2.  It  has  a 
pronounced  maximum  for  j-to/ksT  ~  0.6.  Parameter  nC 3  also  displays  a  maximum,  which  occurs 
for  /acAbT  ~  1.1.  In  contrast,  parameters  nCi,4  depend  on  /Jo/kaT  monotonically. 

The  results  of  Appendices  C  and  D  show  that,  for  the  Lam  A  and  extension  modes  in  crystals  cut 
in  (100)  direction,  the  relative  shift  Scov/(ov  is  determined  by  coefficient  11C2,  which  explains  the 
nonmonotonicity  of  the  shift.  For  crystals  cut  in  (110),  the  shift  of  the  Lame  mode  is  fully 
determined  by  coefficient  nC4  and  is  mono-  tonic,  whereas  for  the  extension  mode  the  expression 
for  the  shift  has  contributions  from  nCi,  nCi,  and  11C4  that  partly  compensate  each  other,  leading 
to  a  comparatively  small  shift  all  together  and  its  weak  nonmonotonicity. 

IV.  Conclusions 

The  results  of  this  paper  show  that  the  electron-  phonon  coupling  strongly  affects  the 
nonlinearity  of  vibrational  modes  in  semiconductor-based  nano-  and  micromechanical  systems. 
The  mechanism  of  the  effect  is  the  strain-induced  redistribution  of  the  electrons  between  the 
valleys  of  the  conduction  band.  The  redistribution  results  from  lifting  the  degeneracy  of  the 
electron  energy  spectrum  by  the  strain  from  a  vibrational  mode.  The  analysis  refers  to  the  range 
of  temperatures  where  the  rate  of  intervalley  scattering  strongly  exceeds  the  frequencies  of  the 
considered  modes.  In  this  case  the  valley  populations  follow  the  strain  adiabatically. 

The  change  of  the  valley  populations  is  a  strongly  non-  linear  function  of  the  strain  tensor.  The 
respective  expansion  of  the  free  energy  in  the  strain  is  an  expansion  in  the  strain  multiplied  by 
the  ratio  of  the  electron-phonon  coupling  energy  (in  particular,  the  deformation  potential)  to  the 
chemical  potential  of  the  electron  system  or  the  temperature.  This  ratio  is  large,  >  103.  It  is  this 
parameter  that  makes  the  nonlinearity  of  the  vibrational  modes  in  doped  semiconductor 
structures  strong. 

Of  special  interest  in  nano-  and  micromechanical  systems  is  the  amplitude  dependence  of  the 
vibration  frequency.  To  the  leading  order,  it  is  determined  by  the  quartic  terms  in  the  expansion 
of  the  free  energy  in  strain.  These  terms  are  comparatively  large  in  doped  crystals.  We  have 
calculated  the  nonlinear  elasticity  tensor  that  describes  the  electron  contribution  to  the  terms  in 
the  free  energy,  which  are  quartic  in  the  strain.  The  explicit  expressions  for  the  tensor 
components  refer  to  semiconductors  with  the  valleys  on  (100)  axes,  in  particular,  to  silicon.  We 
have  also  found  this  tensor  for  germanium.  In  silicon,  along  with  the  deformation  potential 
coupling,  an  important  role  is  played  by  the  coupling  to  shear  strain.  Such  strain  lifts  the  band 
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degeneracy  at  the  zone  boundary  and  is  effectively  described  by  a  two-phonon  coupling.  We 
show  that  this  coupling  also  leads  to  strong  nonlinearity  of  vibrational  modes. 

The  parameter  of  the  electron  coupling  to  shear  strain  in  silicon  is  not  easy  to  access  in  the 
experiment  [26,  32],  Measurements  of  the  nonlinear  frequency  shift  provide  a  direct  means  for 
determining  this  parameter.  In  particular,  the  nonlinear  frequency  shift  of  the  fundamental  Lame 
mode  in  a  silicon  plate  cut  along  (110)  axes  is  determined  by  this  parameter  only,  except  for 
small  corrections  from  the  nonlinearity  of  the  undoped  crystal. 

We  found  that  the  nonlinear  frequency  shift  strongly  depends  on  the  type  of  a  vibrational  mode 
and  the  crystal  orientation.  We  also  found  that  the  ratio  of  the  frequency  shift  to  the  squared 
vibration  amplitude  can  be  profoundly  nonmonotonic  as  a  function  of  electron  density  and 
temperature.  The  results  provide  an  insight  into  the  experimentally  observed  strong  mode 
nonlinearity  in  doped  crystals  [17].  In  terms  of  applications,  they  enable  choosing  the  appropriate 
range  of  doping  and  the  temperature  regime  to  optimize  the  operation  of  nano-  and 
micromechanical  resonators. 
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Appendix  A:  Expansion  of  the  Free  Energy  in  Terms  of  the  Strain-Induced  Shift  of  the 
Energy  Valleys 

The  major  effect  of  a  strain  on  the  electron  free  energy  comes  from  the  shift  of  the  energy 
valleys.  We  will  assume  that  valley  a  is  shifted  in  energy  by  SEa  and  the  shift  is  small, 

\5Ea\  «  max(foT,  /jo),  where  /jo  is  the  chemical  potential  in  the  absence  of  strain.  We  further 
assume  that  the  vibrations  are  slow  compared  to  the  time  it  takes  the  electron  system  to  come, 
locally,  to  thermal  equilibrium  for  given  values  of  SEa,  i.e.,  the  temperature  and  the  chemical 
potential  are  the  same  in  all  valleys.  Since  for  high  electron  densities  the  thermal  conductivity  is 
high,  the  change  of  the  temperature  compared  to  the  ambient  temperature  can  be  disregarded; 
also,  as  men-  tioned  in  the  main  text,  the  electron  density  n  summed  over  all  valleys  is  constant. 


Expanding  the  electron  free  energy  density  to  the  4th  order  in  the  strain-induced  shifts  8Ea,  we 
find  that,  in  an  /V-valley  semiconductor,  the  change  8Fe  of  the  free  energy  density  is 


SFe 

nkuT 


=A€  +  l-l*  [(3, )2  -  AT]  +  to  -  3A*Sf  +  2(3e )31 
2  ri/2  L  i  6/-1/2  L 

+  7  [®2  -  2A? (Ac)2  4-  (Af)4]  + 


1  F"9 

1  M/2 


24  F, 


1/2 


[•1A?  Af  -  Aj  -  GA2  (AE)2  +  3(Af)4]  . 


(At) 
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Here,  A™  =  AT_1  'y\a(SEa/kyjT)n'1 .  We  use  the  standard  notation 

Fi/2  (z)  =  /o'0  dyy1^1  f\  1  +  exp(y  —  t)]  primes  dFWdx.  Function  Fm  and  its  derivatives 
are  calculated  for  x  =  hoAbT. 

Equation  (Al)  immediately  gives  the  tensors  Anof  the  expansion  of  the  free  energy  increment 
(6)  if  one  expresses  the  shift  5Ea  of  the  valleys  in  terms  of  the  strain  tensor.  In  the  deformation 
potential  approximation  the  relation  between  SEa  and  F  is  given  in  the  main  text,  see  also 
Eq.  (A2)  below. 

In  the  case  of  Si  crystals,  which  are  often  used  in  micromechanical  resonators,  an  important 
contribution  to  SEa  comes  from  the  shear-strain  induced  splitting  of  the  electron  energy  bands  at 
the  zone  boundary.  Shear  strain  does  not  lead  to  the  valley  shift  in  the  deformation  potential 
approximation.  The  overall  shift  of  valley  a,  to  the  lowest  order  in  the  coupling  that  causes  it 
(i.e.,  to  the  first  order  in  the  deformation  potential  where  its  contribution  is  nonzero  and  to  the 
second  order  in  the  band  splitting  for  shear  strain)  is  [26] : 

SEa  =  £  3£'e«  -  3*4,  5*  =  ^-.  (A2) 

ij 


Here  we  use  that  silicon  has  six  valleys  located  at  the  (100)  axes,  and  we  chose  the  coordinate 
axes  x,  y,  z  along  (100).  Respectively,  the  valley  index  a  takes  on  three  values  that  correspond 
to  the  x,  y,  z  axes  (the  valleys  lying  on  the  same  axis,  but  in  the  opposite  directions,  are 
equivalent).  The  strain  sa,  which  enters  the  second  term  in  the  right-hand  side  of  Eq.  (A2),  is  a 
component  of  the  strain  tensor  s,j  with  i,  j  such  that  i,  j  F  a  and  i  F  j.  The  parameter  2 EU’  is  the 
interband  matrix  element  of  the  electron-phonon  coupling  calculated  for  the  electron  conduction 
bands  Ai  and  A2’  at  the  X  point  on  the  boundary  of  the  Brillouin  zone,  where  the  bands  cross;  A E 
is  the  energy  separation  between  the  bands  Ai  and  A2’  at  the  value  of  the  wave  vector  k  that 
corresponds  to  the  conduction  band  minimum.  Parameter  ESh  is  the  effective  deformation 
potential  of  two-phonon  coupling  to  shear  strain.  The  numerical  value  of  Ssh  is  not  well  known. 

The  experimental  data  give  S  u’  ~  7  -  8  eV  [26,  32]  and  the  numerical  data  on  the  band  splitting 
give  A E  ~  0. 7  eV  [33]  so  that  Ssh  is  in  the  range  of  280  -  360  eV;  this  is  essentially  an  order  of 
magnitude  estimate. 

In  calculating  SFe  in  Eq.  (Al)  we  kept  terms  that  are  quartic  in  F .  The  components  of  the  tensors  A& 
in  Eq.  (6)  are  expressed  in  terms  of  SFe  as 


SkSF, 


OSj1jl  ...OEikjk 


(A3) 


Tensors  A  are  symmetric  with  respect  to  the  interchange  of  indices  4  jk  and  the  pairs 
iikjk)  <->■  0V  jk).  For  the  considered  long-wavelength  strain,  tensors  A;,  are  independent  of 
coordinates.  The  corrections  A2  to  the  linear  elasticity  tensors  were  found  previously  [19,  20] 
and  are  not  discussed  in  this  paper. 
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Appendix  B:  Nonlinear  Elastic  Constants  of  Germanium 


In  this  section  we  provide  the  corrections  to  the  nonlinear  elastic  constants  of  germanium,  which 
are  due  to  the  redistribution  of  the  electrons  over  the  valleys.  Germanium  has  four  equivalent 
valleys  in  the  conduction  band,  which  are  located  on  the  boundary  of  the  Brillouin  zone  along 
(111)  axes.  We  use  the  Voigt  notation  and  write  the  components  of  the  corrections  to  the 
nonlinear  elasticity  tensor  Sc ~  in  the  frame  where  the  axes  (x,  y,  z)  are  along  the  (100)  directions 
of  the  crystal.  Using  the  results  of  Appendix  A,  we  obtain 
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<Sc4444 
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T7-3  F" 

”  ii  "  1/2 
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The  notations  are  the  same  as  in  Appendix  A  and  in  Table  I.  The  electron-phonon  coupling  does 
not  contribute  to  the  other  third-  and  fourth-order  elastic  constants. 

Corrections  SC44  and  SC456  for  germanium  were  found  by  Keyes  [18];  however,  his  final 
expression  for  Sc456  differs  from  Eq.  (Bl)  by  a  factor  of  4  (our  expressions  for  Sc44  coincide  with 
Ref.  [18]).  Parameters  JC4444  and  <5c4455  have  not  been  found  before,  to  the  best  of  our  knowledge. 
In  the  limiting  cases,  corrections  dc4444  and  Sc4455  have  the  same  dependence  on  temperature  and 
electron  density  as  constant  nCi  discussed  in  Sec.  IV  A. 

Appendix  C:  Duffing  Nonlinearity  Parameter  for  a  Lame  Mode  in  a  Square  Single- 
Crystal  Plate 

We  consider  a  square  plate  with  side  L  and  thickness  h  made  out  of  a  single  crystal  with  cubic 
symmetry.  If  the  crystal  is  cut  out  along  (100)  or  (110)  axes,  one  of  the  simplest  modes  is  the 
first  Lame  mode  [34],  The  normalized  displacement  field  is 

=  \/2cos(7raT/L)sin(7ry/L), 
u[y'!  =  —  \/2sin(7r2’/L)  cos(7r  y/L).  (Cl) 


Here,  x  and  y  axes  are  in  the  lateral  plane  along  the  sides  of  the  square,  axis  z  is  perpendicular  to 
the  plate  and  =  0,  Calculating  the  strain  tensor  for  the  displacement  (Cl)  and  substituting 
the  expressions  into  Eqs.  (8)  and  the  relation 

Muil  =  f  *A<f>  •  i(i,)  ®  ew\  (C2) 
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for  the  plate  cut  out  along  (100)  axes  we  obtain,  in  Voigt  notation  for  the  elasticity  tensors, 


Mul  =  7 r2/t  (cu  -C12), 


lv  = 


3tt4/? 

igl2 


(cun  —  4C1112  +  3^1122)  • 


(C3) 


If  we  consider  silicon  and  take  into  account  only  the  contribution  Sc  to  the  nonlinear  elasticity 
tensor  c",  with  the  account  taken  of  Table  I,  the  expression  for  yv  simplifies  to 

7„  =  (27w4h/32L3)6cnu.  (C4) 

For  the  Lame  mode  cut  along  the  (110)  axis,  if  the  tensors  are  calculated  in  the  axes  (100),  we 
have 


1 =  2n2hc44, 

lv  =  (37r4/i/2L2)<5c4444.  (C5)i 

Note  that  only  coupling  to  shear  strain  contributes  to  the  nonlinearity  parameter  yv  in  this  case. 

Appendix  D:  Duffing  Nonlinearity  Parameter  for  an  Extension  Mode  in  a  Single-Crystal 
Narrow  Beam 


We  consider  the  fundamental  extension  mode  in  a  thin  beam  of  length  L  with  a  rectangular  cross- 
section  of  area  S  «  L2.  The  beam  is  cut  along  a  symmetry  axis,  and  the  sides  are  also  along 
symmetry  planes  of  a  cubic  crystal.  From  the  free-surface  boundary  conditions,  the  normalized 
displacement  field  is  [34]: 


u 


u 


y 

(*d 


\/2cos(7t  x/L), 

y/2n<r2  .  ,  /r, 

—  - ysin(-Kx/L) 

.  ,  /r, 

—  - 2  sin(7Tj:/L) 


(Dl) 


This  expression  takes  into  account  transverse  compression  that  accompanies  beam  extension  and 
uses  the  smallness  of  the  beam  cross-section;  corrections  ~  S/L 2  are  disregarded.  The  transverse 
compression  in  a  cubic  crystal  cut  in  a  symmetric  direction  is  described  by  Poisson’s  ratios  02 
and  (73.  Generally,  they  do  not  coincide.  In  Eq.  (Dl)  the  transverse  coordinates  y  and  z  are 
counted  off  from  the  center  of  the  beam 


For  the  longitudinal  direction  of  the  beam  (100)  and  the  sides  parallel  to  (100)  planes,  the 
Poisson  parameters  are  equal,  02  =  03  and  a  =  02  =  <73  =  c\t/(c\\  +  cn).  In  this  case  Eqs.  (8)  and 
(C2)  give 
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n2S  (cn  (cn  +  c12)  -2c?2) 


L  (cn  +  c12) 


=  (7r45/4L3)  [ci in  —  8o-cih2 
+  12<T2(ch22  +  Cn23)  —  8<T3(Cm2  +  3cn23) 
+2cr4(cmi  +  4cm2  +  3ch22)]  .  (D2) 


The  expression  for  yv  is  simplified  if  in  the  nonlinear  elasticity  tensors  we  take  into  account  only 
the  contribution  from  the  electron-phonon  coupling  as  given  in  Table  I  and  also  allow  for  the 
interrelation  between  different  components  of  the  tensor  8c".  Then  for  a  silicon  beam 


lv  =  (tt45/4L3)(1  +  o-)4d'cuu.  (D3) 


For  extension  along  (110)  axis,  with  one  side  parallel  to  (100)  plane  and  the  other  side  parallel  to 
(1 10)  plane,  the  Poisson’s  ratios  02  =  cr(l  10,  110)  and  02  =  cr(l  10,  001)  are  given  in  Ref.  [35]. 
Then  Eqs.  (8)  and  (C2)  give 

^  2  _  47t2S  C44  (cii(cn  +C12)  —  2c?2) 

L  cn(cn  -I-  ei2  +  2044)  -  2 cf2 
7r*S  r 

lv  =  22]FT  pill  (<*2  -  4o2  +  G<t2  -  4<t2  +  8*i  +  1) 

+  4cni2  (<T2  —  1 )  (<7o  +  2<7^<73  -  3<T2  -  4(72^3  +  3(72  +  803  +  2(73  -  l) 

+  3cii22  (<72  —  l)2  (<72  —  2(72  +  8(73  +  l)  +  24<?1123<73  (<72  —  l)2  (<72  +  <73  —  1) 

-I-  48cil44<73  (<72  +  l)2  +  9Gcj244<73  (<72  -  1)  (<72  +  l)2  +  24cil55  (<72  —  l)2 


+24CJ2CC  (<72  —  l)  +  8C4444  (<72  +  l)4j  . 


(D4) 


If  in  the  nonlinear  elasticity  tensor  c"  we  take  into  account  only  the  contribution  8c"  from  the 
electron-phonon  coupling,  in  the  case  of  a  silicon  beam  the  expression  for  yv  simplifies  to 


Expressions  (D2)  and  (D4)  were  generated  using  a  computer  code  to  calculate  the  sums  and 
integrals  in  Eq.  (8). 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 


ACRONYM 

AFRL 

DARPA 

Ge 

MEMS 

Si 


DESCRIPTION 

Air  Force  Research  Laboratory 
Defense  Advanced  Research  Agency 
Germanium 

Microelectromechanical  Systems 
Silicon 
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